function [x, y, a, blink, speed, vel_x, vel_y, shrink]=EL1k_blinkremover_2013(t, x, y, a, thresh, plot_it)
%[x1 y1 a1 blink speed vel_x vel_y shrink]=EL1k_blinkremover_2013(t, x, y, a, [thresh], [plot_it])
% data must correct to degrees per pixel first.
% inputs:
% t = time in ms
% x = horizontal position in degrees of visual angle
% y = vertical position in degrees of visual angle
% a = area : size of pupil
% thresh = threshold for cutting out 'blink' data and replacing it with interpolated lines
% plot_it = do you want to see my superfluous ploting that takes eons?
% ~bcoe:  coe@queensu.ca

% 2013: updated start loc and end lo arrays and removed 2 the threshold catch
% also added the ability to skip short bits os data (sss-eloc(i-1)<5)
%
if nargin==0
    help(mfilename);
    return
end
if nargin<4
    warning ([mfilename ' requires 4 input vectors (t, x, y, a) of equal length']) %#ok<WNTAG>
    help(mfilename);
    return
end
if nargin==4
    thresh=2.5;
    plot_it=0;
end
if nargin==5
    plot_it=0;
end

vel_x=single(zeros(size(x)));
vel_y=single(zeros(size(x)));
speed=single(zeros(size(x)));
shrink=single(zeros(size(x)));
%blink=false(size(x));
if sum(isnan(x))==length(x) || sum(isnan(y))==length(y)
    blink=false(size(x));
    return
    %             vel_x()=0;
    % 			vel_y()=0;
    % 			speed(starloc(tr):end_loc(tr))=0;
    % 			shrink(starloc(tr):end_loc(tr))=0;
    % 			x(starloc(tr):end_loc(tr))=x;
    % 			y(starloc(tr):end_loc(tr))=y;
    % 			a(starloc(tr):end_loc(tr))=a;
    % 			blink(starloc(tr):end_loc(tr))=1;
    % 			return
end
blink=a<(mean(a(a>0))/2); % less than half of it's non-zero mean
if isnan(sum(x))
    tdiff=[(t(2)-t(1));(t(3:end)-t(1:end-2))/2;(t(end)-t(end-1))];
    adiff=[(a(2)-a(1));(a(3:end)-a(1:end-2))/2;(a(end)-a(end-1))];
    tshrink=single(filtfilt([1 1 1]/3,1,double(abs(adiff./tdiff)))); % filtfilt now requires doubles
    tshrink(isnan(x))=nan;
    [val,cnt]=int_hist(tshrink);ccnt=cumsum(cnt)/sum(cnt);
    %val(find(ccnt>.90,1,'first'))
    tshrinkstd=std(tshrink(tshrink<val(find(ccnt>.95,1,'first'))));% temp var just to shore up the ends
    firstgood=find(tshrink<tshrinkstd  &~blink & ~isnan(x) ,1,'first');
    lastgood=find(tshrink<tshrinkstd  &~blink & ~isnan(x) ,1,'last');
    if firstgood>1
        x(001:firstgood)=x(firstgood);
        y(001:firstgood)=y(firstgood);
        a(001:firstgood)=a(firstgood);% why was this zero
        %tshrink(001:firstgood)=0;
    end
    if lastgood<length(x)
        x(lastgood:end)=x(lastgood);
        y(lastgood:end)=y(lastgood);
        a(lastgood:end)=a(lastgood);% why was this zero
        %tshrink(lastgood:end)=0;
    end
    %a(1)=a(find(isnan(a)==0,1,'first'));
    %a(end)=a(find(isnan(a)==0,1,'last'));
end %isNaN
blink=a<(mean(a(a>0))/2); % less than half of it's non-zero mean

[speed, vel_x, vel_y, shrink]=xyt2speed(x,y,t,a,3);

if plot_it
    aNAN=a;aNAN(a==0)=NaN;%just for plotting
    aNAN=(aNAN-nanmean(aNAN));
    aNAN=(aNAN/nanstd(aNAN))*5;% mean of zero std of 10;
    scaleit=10;
    ssr=mean(speed(speed<nanmean(speed)*2))/mean(shrink(shrink<nanmean(shrink)*2));% speed:shrink ratio
    ttt=t-t(1); %just for plotting
    figure(1);clf; hold on;
    ylim([-20 20])
    xlim([0 ttt(end)])
    hh=plot(ttt,aNAN,'color',[.5 .5 .5],'linewidth',2);
    hh(2)=plot(ttt,speed/scaleit,'.-g');
    hh(3)=plot(ttt,shrink/scaleit*ssr,'.-','color',[1 .5 0]);
    hh(4)=plot(ttt,x,'b','linewidth',2);
    hh(5)=plot(ttt,y,'r','linewidth',2);
    
end

s_blink01=filtfilt([1 1 1 1 1 1 1]/7,1,double(blink));% filtfilt now requires doubles
s_blink02=filtfilt([1 1 1 1 1 1 1]/7,1,round(s_blink01+.25));% filtfilt now requires doubles
blink=logical(round(s_blink02+.25));
blink=blink | (isnan(x) | isnan(y)); % 

blink(end)=0;
%loc=find(diff([0; blink]))-1;
loc=(diff([0; blink]));
loce=[find(loc==1) find(loc==-1)-1];

spinner={'\','|','/','-'};
fprintf('running blink remover:-\n');
if ~isempty(loc) %actually blinked
    shrinkstd=std(shrink(shrink<100))*thresh;% was 10;
    speedstd=nanstd(speed(speed<200 & speed>0))*thresh;
    speedstd=inf;
    
    if plot_it
        plot(ttt(blink),aNAN(blink),'xk');
        vfill(ttt(loce(:,1)),ttt(loce(:,2)),[],.2 )
        %vfill(ttt(loc(1:2:end)),ttt(loc(2:2:end)),[],.2 )
        vline(ttt(loce))
        %hline(shrinkstd/5,'c');
        hline(speedstd/scaleit,'-.g');
        hline(shrinkstd/scaleit*ssr,[1 .5 0]);
    end
    
    sloc=[1; loce(:,1);length(blink)];
    eloc=[1; loce(:,2);length(blink)];
    for i = 2:length(sloc)-1
        fprintf('\b\b%s\n',spinner{mod(i,4)+1});
        %forward in time
        eee=find(speed(eloc(i):sloc(i+1))<speedstd &...
            abs(shrink(eloc(i):sloc(i+1)))<shrinkstd ...
            , 1,'first')+eloc(i)-1 ;
        if isempty(eee)% % try with next blink
            eee=eloc(i+1);
            sloc(i+1)=sloc(i);
        end
        
        if plot_it && eee
            xlim([ttt(eloc(i))-50 ttt(eloc(i))+150 ])
            vfill(ttt(eloc(i)), ttt(eee),[],.5)
            %drawnow
        end
        eloc(i)=eee;% redefine end
        
        % backward in time
        sss=find(speed(eloc(i-1):sloc(i))<speedstd & ...
            abs(shrink(eloc(i-1):sloc(i)))<shrinkstd ...
            ,1,'last')+eloc(i-1)-1;
        if isempty(sss)% use previous blink
            sss=sloc(i-1);
            %sloc(i+1)=sloc(i-1);
        end
        if sss-eloc(i-1)<5 && sss~=sloc(i-1); % if gap of real data too small
            drawnow
            sss=sloc(i-1);
        end
        if plot_it
            xlim([ttt(sloc(i))-150 ttt(sloc(i))+50 ])
            vfill(ttt(sss), ttt(sloc(i)),[],.5)
            %drawnow
        end
        sloc(i)=sss;% redefine start
        
        % interpolate
        %try
        %         if length(sss:eee)>500
        %             drawnow
        %         end
        blink(sss:eee)=1;
        x(sss:eee)=linspace(x(sss),x(eee),eee-sss+1);
        y(sss:eee)=linspace(y(sss),y(eee),eee-sss+1);
        a(sss:eee)=linspace(a(sss),a(eee),eee-sss+1);
        % recalculate
        [speed, vel_x, vel_y, shrink]=xyt2speed(x,y,t,a,3);
        if any(isnan(x(1:eee))| isnan(y(1:eee)))
            drawnow
            badapple=find(isnan(x(1:eee)));
            if isempty(diff(badapple))% only 1
                x(badapple)=mean(x([badapple-1,badapple+1]));
                y(badapple)=mean(y([badapple-1,badapple+1]));
                a(badapple)=mean(a([badapple-1,badapple+1]));
            else
                if all(diff(badapple)==1)
                    x(badapple)=linspace(x(min(badapple)-1),x(max(badapple)+1),range(badapple)+1);
                    y(badapple)=linspace(y(min(badapple)-1),y(max(badapple)+1),range(badapple)+1);
                    a(badapple)=linspace(a(min(badapple)-1),a(max(badapple)+1),range(badapple)+1);
                else
                    warning('bad apples in blinkremover!')
                end
            end % if isempty(diff(badapple))
        end% if any(isnan(x(1:eee))| isnan(y(1:eee)))

        if plot_it
            aNAN(sss:eee)=linspace(aNAN(sss),aNAN(eee),eee-sss+1);
            % replot
            xlim([ttt(sss)-150 ttt(eee)+150 ] )
            %my_range=
            
            set(hh(1),'ydata',aNAN)
            set(hh(2),'ydata',speed/scaleit)
            set(hh(3),'ydata',shrink/scaleit*ssr)
            set(hh(4),'ydata',x)
            set(hh(5),'ydata',y)
            
            plot(ttt(sss:eee), x(sss:eee),'ob')
            plot(ttt(sss:eee), y(sss:eee),'or')
            plot(ttt(sss:eee), aNAN(sss:eee),'ok')
            plot(ttt(sss:eee), speed(sss:eee)/scaleit,'og')
            plot(ttt(sss:eee), shrink(sss:eee)/scaleit*ssr,'o','color',[1 .5 0])
            drawnow
            pause (.1)
        end
        %catch %#ok<CTCH>
        %	drawnow;
        %end
    end % sloc
end % actually blinked

end % function



